#ifndef E3D_MATRIX4_HPP
#define E3D_MATRIX4_HPP

#ifndef __cplusplus
#error "C++ required"
#endif

namespace e3d {

	template <typename real>
	class Matrix4 {
	public:
		typedef real_traits<real> traits;
		typedef Matrix4<real> matrix_type;
		typedef Vector3<real> vector_type;

	private:
		real v[16];
	public:
		Matrix4() : v({
			1, 0, 0, 0,
			0, 1, 0, 0,
			0, 0, 1, 0,
			0, 0, 0, 1}) {}

		Matrix4(Matrix4 const & o) {
			for (int i = 0; i < 16; ++i)
				v[i] = o.v[i];
		}

		matrix_type & operator *= (const matrix_type & o) {
			return *this;
		}

		real const & operator[](int index) const {
			if (index < 0 || index >= 16)
				throw outofbound("matrix access out of bound");
			return v[index];
		}

		matrix_type & transpose();
		matrix_type & invert();

		static matrix_type transposed(matrix_type const & m) {
			matrix_type r(m);
			r.transpose();
			return r;
		}

		static matrix_type inverted(matrix_type const & m) {
			matrix_type r(m);
			r.invert();
			return r;
		}

		static matrix_type translation(vector_type const & v) {
			matrix_type m;
			m.v[3] = v.x;
			m.v[7] = v.y;
			m.v[11] = v.z;
			return m;
		}

		static matrix_type rotation(vector_type const & v) {
		}

		static matrix_type scaling(vector_type const & v) {
			matrix_type m;
			m.v[0] = v.x;
			m.v[5] = v.y;
			m.v[10] = v.z;
			return m;
		}

		template <typename R>
		friend 	Vector3<R> operator*(Matrix4<R> const &, Vector3<R> const &);

	};

	template <typename R>
	Vector3<R> operator*(Matrix4<R> const & m, Vector3<R> const & v) {
		R a, b, c, d;
		a = m[0] * v.x + m[1] * v.y + m[2] * v.z + m[3];
		b = m[4] * v.x + m[5] * v.y + m[6] * v.z + m[7];
		c = m[8] * v.x + m[9] * v.y + m[10] * v.z + m[11];
		d = m[12] * v.x + m[13] * v.y + m[14] * v.z + m[15];
		return Vector3<R>(a / d, b / d, c / d);
	}

}

#endif
